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CN , We give a full account of the Numerical Stochastic Perturbation Theory 

^ | method for Lattice Gauge Theories. Particular relevance is given to the in- 

clusion of dynamical fermions, which turns out to be surprisingly cheap in 
this context. We analyse the underlying stochastic process and discuss the 
convergence properties. We perform some benchmark calculations and - as a 
byproduct - we present original results for Wilson loops and the 3-loop critical 
^"T 1 ■ mass for Wilson fermions. 
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1. Introduction 

Perturbation Theory applied to Quantum Field Theories in Lattice regularization is very 
difficult. One has to face the difficulty of handling complicated trigonometric functions 
defined in the Brillouin zone (the momentum space on the lattice). Things become even 
more cumbersome in the case of Lattice Gauge Theories due to appearance of new Feynman 
vertices at any new perturbative order. This makes computations extremely difficult at 2- 
loops level and virtually unfeasible beyond that. Ideally Perturbation Theory would not be 
necessary on the lattice. However, in practice, it is needed for a number of reasons: to match 
perturbative results obtained in a continuum regularization (typically MS) with the non- 
perturbative ones from the lattice; to compute perturbative renormalization factors of bare 
parameters and operators; to determine the so-called improvement coefficients for lattice 
actions and operators. In recent years much progress has come from a non-perturbative 
implementation of many of these tasks. Still, perturbative results are needed when non- 
perturbative results are not available (and difficult to achieve). Also, one would like to 
gain insight by comparing perturbative and non-perturbative results (and a fair comparison 
actually calls for going beyond one loop). In the end, one knows that Perturbation Theory 
has to be reliable in a given limit of the theory and it is important to carefully assess what 
the results are in this limit. For a recent review of Perturbation Theory on the lattice see 
PP. A pioneering high order (three loops) computation in QCD with Wilson fermions can 
be found in j2|. 

Numerical Stochastic Perturbation Theory (NSPT) was introduced E] as a numerical 
application of Stochastic Quantization E] and successfully applied [3 El El HHj to perform 
high order perturbative calculations in Lattice Gauge Theories (LGT). The basic idea is to 
integrate on a computer the differential equations of Stochastic Perturbation Theory. This 
results in a slight modification of a non perturbative Langevin algorithm. Until very recently 
the applications of NSPT to LGT were limited to the quenched case. Given the similarity 
between NSPT and a non-perturbative Monte Carlo, one could fear that unquenched NSPT 
would cost some order of magnitude more than the quenched case, just like the non pertur- 
bative algorithm from which it is derived. The main purpose of this paper is to describe 
the unquenched version of NSPT for LGT and present some benchmark analysis. The good 
message is that full NSPT is not much more expensive than quenched NSPT, provided a 
fair efficiency is granted for a basic tool, namely the Fast Fourier Transform. In fact the 
perturbative expansion of the inverse fermion matrix only requires the inverse of the free 
fermion matrix. We will show how to efficiently exploit this property. 

Since the first introduction of NSPT in 011], much experience has been collected which 
brought to a better understanding of the underlying stochastic process. The interest about 
the stochastic process was strongly motivated by the observation that, when applied to very 
simple (0-dim) models, NSPT led to very large (and not normally distributed) fluctuations 
at high orders. These models were studied in [TT]. Here we consider specifically the case 
of LGT. In order to rely on NSPT as a numerical method to extract physical results, one 
needs to know under which conditions the stochastic process converges to a limit distribution 
and assess - as much as possible - the properties of the limit distribution and the rate of 
convergence. We will see that a limit distribution exists, if some kind of gauge fixing has been 
introduced and if the fields do not contain zero modes. The assessment of the properties of 
the limit distribution - which is crucial - was already discussed in [TTJ. There it was shown 
that LGT do not display the pathological behaviour of the O-dimensional models, and an 
argument was given to understand why this happens. As a benchmark for error analysis, in 
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this paper we collect results for a large set of Wilson loops and show how the fluctuations 
depend on the loop size and the perturbative order. 

It is often useful to compute gauge dependent quantities in a fixed gauge. In the context 
of NSPT a form of gauge fixing is always imposed by a stochastic term inspired to the gauge 
fixing procedure introduced by Zwanziger in ^2] (actually our prescription is a perturbative 
form of [HD- However, it is not clear what is the perturbative relation between this gauge 
fixing and the more popular ones. Still, we can show a simple way to recover Landau gauge 
fixing. We can also perform computations in covariant gauges without introducing ghost 
fields. The Faddeev-Popov determinant is represented by the same technique used for the 
fermionic one. 

Another interesting feature of NSPT is the possibility of computing, with the same effort, 
perturbative expansion around non trivial vacua. The vacua configurations may even be 
known only numerically. 

In section 2 we briefly review how NSPT can be naturally derived from the context of 
Stochastic Quantization and how this can be specified to the case of LGT. In section 3 we 
concentrate on the analysis of the NSPT stochastic process for LGT: we first explain why a 
stochastic gauge fixing is needed and how it can be practically implemented; then we show 
how the modes with zero momenta must be regularized in this context. These corrections 
allow us to prove the convergence of the stochastic process. We conclude the section with a 
comment on simulation times and their dependence on the lattice volume and the maximum 
perturbative order. In the same spirit of assessing the efficiency of the method we also 
perform the error analysis of a large set of Wilson loops in the quenched approximation, 
discussing the dependence on the loop size and the perturbative order. Section 4 is devoted 
to the illustration of some special applications. First we show how to evaluate in perturbation 
theory gauge dependent quantities in Landau gauge. We also comment on the strategy for 
the implementation of the Faddeev-Popov mechanism for general covariant gauges. Some 
more comments on this subject are differed to the following discussion of fermions, since it 
relies on the same technique. We then consider perturbative expansions around non trivial 
vacua. Section 5 is devoted to the inclusion of dynamical fermions, while in section 6 we 
present benchmark unquenched computations for rif = 2 Wilson fermions, i.e. Wilson loops 
and the critical mass to the third order. In each section we also provide directions for a 
practical implementation of the method in a Monte Carlo simulation. In section 7 we give 
our conclusions. 

2. From Stochastic quantization to NSPT for Lattice Gauge Theories 

2.1. Stochastic Quantization. NSPT was introduced [3J H] as a numerical application of 
Stochastic Quantization It is well known that from the latter one can formulate a 

Stochastic Perturbation Theory, which converges (in a sense that will be clear in a moment) 
to the standard perturbative expansion of field theories. NSPT is nothing but the numerical 
implementation of this program. In order to fix the main ideas and introduce the notation let 
us remind the fundamental assertion of Stochastic Quantization. One starts with an action 
for a field theory S[4>] and aims at computing expectation values with respect to the path 
integral measure, i.e. 



The basic fields of the theory are now given an extra degree of freedom t {4>{x) h- > tfi^x^t)). 
One can think of it as a fictitious (stochastic) time in which an evolution takes place according 
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to the Langevin equation 

d^jx-t) _ dS[4>] 

[ ) dt ~ dM^t) v[ , ) ' 

The first term is the drift term given by the equations of motion (which would drive the field 
to the classical solution). One adds to the latter a Gaussian noise 

rj(x; t) : (r](x, t) rj(x', t')) v = 2 5(x- x') 5(t - t') 

where 

f Dr](z,T) ... e -V dzdr V 2 (z,r) 

" ~ j D V {z,T)e-y dzdT ^) 

77 is the source of the stochastic nature of the process (this motivates the notation (f) v ). The 
natural expectation values to compute in this framework are with respect to this noise. The 
main assertion of Stochastic Quantization is that [3] 

(2) (0[^(x x ; t) ... (f) v {x n \ t)])^^^oo(O[0(xi) . . . (f>{x n )]), 

which means that in the limit of the stochastic time going to infinity expectation values of 
the field with respect to the Gaussian noise converge to the functional integral expectation 
values one is interested in 1 . One can now move to Stochastic Perturbation Theory (SPT) 
by noting that S[<f>] can be written as a sum S[4>] = Sq[</>] + Sj[<p,g}, in which an interaction 
part Sj (which is a function of a coupling constant g) has been singled out. For the free 
field action So[<f>] an integral solution for the Langevin equation exists, which is written 
(often in momentum space) in terms of a Green function. For the full theory, the Langevin 
equation can be converted into an integral equation, whose iterative solution as a series in the 
coupling constant actually results in SPT. This is the original approach to SPT, for which 
the interested reader can refer to jHJ. Notice that this same approach can be recovered by 
looking for a solution of the Langevin equation as a (formal) series in the coupling constant 

(3) J? (x;t) = ^(x;t) + ^^(x;t). 

n>0 

If one plugs this expansion in the Langevin equation, the latter gets translated into a hier- 
archy of equations, which one can exactly truncate at any given order. Their solutions are 
the same integral expressions mentioned above. In the simple case of a bosonic non-gauge 
theory 2 it is possible to show that stochastic diagrams are generated, which can be evaluated 
in the t — > 00 limit. In this limit they reconstruct the contribution of the standard Feynman 
diagrams |14j . Another, more formal, proof of the equivalence of SPT and standard field 
theoretic Perturbation Theory relies on the formalism of the Fokker-Planck equation 3 : we 
will sketch very briefly what is the content of to which we refer the interested reader. Let 
us first of all remind that the fundamental assertion contained in Eq. (J2J) can be rephrased in 
terms of a distribution function of the field. The first step is the introduction of a stochastic 
time dependent distribution function for the field according to 

1 Notice that all the fields are evaluated at the same stochastic time in Eq. 
2 Notice that diagrammatic arguments are not conclusive in the case of gauge theories. 
3 Within this framework the gauge theories case is neatly treated once the Stochastic Gauge Fixing term 
has been added. 
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The Langevin equation Eq. can then be traded for the so-called Fokker-Plank equation 
which expresses the (stochastic) time derivative of P[4>,t] 

The essential steps taken in ^2] can be summarized as in the following. 

• P[4>,t] is expanded as a power series in the coupling of the theory 

P[<M] = ]T(M[<M]. 

fc=0 

This also turns the Fokker-Plank equation into a hierarchy of equations. By inspect- 
ing these, the following results can be obtained. 

• One can prove that Po[4>, t]— >t^ooPo q [4>] = 6 , i-e. Po[(f),t] converges to the free 
field path integral measure. 

• In a convenient weak sense it is also true that P k [(f>,t]^ t ^ 00 P^ q [(/)]. 

• The various P k 9 [4>] are related by a set of relations in which one recognizes the 
Schwinger-Dyson equations. Since the solutions of the Schwinger-Dyson equations 
are unique in Perturbation Theory, one has recovered the standard field theoretic 
perturbative expansion. 

• In the case of gauge theories, a key role in getting through the steps sketched above 
is played by the so called Stochastic Gauge Fixing. 

This ends our introduction to Stochastic Quantization and Stochastic Perturbation The- 
ory. As already said, NSPT is a numerical implementation of SPT, which in particular we 
apply to LGT. The numerical nature opens the door to all the problems of numerical sta- 
bility. It also motivates a peculiar treatment of the question of convergence properties of 
NSPT stochastic process for LGT, which we will address in section 3. As motivated in the 
introduction, this issue is a practically important one. This will also lead us to the intro- 
duction of Stochastic Gauge Fixing, which, as said, is crucial for the results of [TH] . We now 
show the practical implementation of NSPT for LGT. 



2.2. NSPT for Lattice Gauge Theories. Consider the Euclidean Wilson action (for the 
gauge group SU(N C )): 



8 ° = -wXH u * +u - 



Just as before, Stochastic Quantization amounts to considering a set of fields U x ^(t; rf) which, 
besides the usual dependence on space-time (x G M 4 ) and direction (/x = . . . 3), also depend 
on the stochastic time i el and (parametrically) on a random field rj. The Langevin equation 
for the U fields now reads 

(6) ti Ux ^ t] v) = Hv* M s G [£/] - Mfc M (t)) £M*; v), 
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where V XAt = T a V a x/1 = T a V 

u xu 1S a l e ft derivative on the group. We recall the definition of 
the Lie derivative (V is an element of the relevant Lie group) 

V a v f{V) = \\m -{f{d aTa V) - f(V)), 

whose relevant property is that integration by part is admitted. The T a are the hermitian 
generators of the algebra ([T a ,T b ] = if abc T c and Tr(T a T b ) = \5 ab ). The T](t) = T a i] a (t) are 
now random fields with a Gaussian distribution that satisfies 4 

(V a (t)) v = (t&(t)i&,(0>„ = 2 5 ab 5„ IM 8 yx 5(t-t'), 

and higher cumulants vanish. Again, one can prove El that under these assumptions the 
gauge fields distribute according to the measure e~ s ^ in the limit of large t. In other words 
(if U(t; if) is a solution of © determined by if): 



]im{0[U(t',rj)])„ = \\ DUe~ SG ^0[U]. 

t^oo 1 Z J 



Having in mind a numerical integration on a computer, the Langevin time can be dis- 
cretized (with step e and t = ne). Following [TH], one can check that a solution to © is 
given by 

(7) U Xft (n + 1; v) = e~ F '^ U m (n; V ) 
where 

(8) I'xft 

[U,rj\ = eV xlx S G [U] + Vei] x/1 



(9) 

and the matrix degrees of freedom of the random field rj are subjected to 5 



(Vi,k( Z ) Vl,m( W )) v 



&il ^ km $ik $lrn 
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Eq. (0) is basically an Euler scheme for Eq. (jUJ) taking care of not leaving the group manifold. 
Being an Euler scheme, results are to be extrapolated linearly as e — > in order to recover 
a solution of Eq. (0). Notice also that the drift F XfJL is a local expression: in order to update 
the link U X}X it suffices to compute the plaquettes Up insisting on U Xil . 

One can now proceed as stated above. Since VSq obviously depends on the coupling, also 
the fields {U Xfl (t; rf)} acquire a dependence on j3 through the Langevin equation ©. As a 
consequence we can write, at least formally, 

(io) m^-i+E^ 72 ^?^)- 

k=l 

The expansion starts with the unity operator. This is actually a choice for the vacuum 
configuration around which we are going to compute perturbative corrections. As we will 
see later, this is not the only possible choice. Due to the presence of a (3 in front of the 
Wilson action, if one now plugs the expansion (|1U|) in F, the latter starts at order yf]3 for 



4 We will often omit the space-time and direction indices as well as explicit rj dependence whenever it is 
reasonable to assume that no confusion can be made. For instance U(t) = U x ^{t\ rj), and r] a (t) = rj^Jt). 

5 In this formula the focus is on matrix indices, z and w are multi-indices collecting the space and 
(stochastic) time degrees of freedom. 
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the drift term, while has a zero order contribution from the random field rj. This makes 
a perturbative evaluation of Eq. (JJJ) inconsistent. However we can redefine the time step 
e' = e/3, in such a way that the first non zero contribution to F is of order ^ and it comes 
both from the drift and from the random noise rj. 

The introduction of the expansion (jl(J|) inside the Langevin equation transforms (J7|) into 
a system of equations, one for each perturbative component of the field 6 : 

(11) UW' = UW-F® 

Jj(2)' = Jj(2) _ F (2) _j_ l_p(l)2 _ ^(1)^(1) 

£/(3)' _ [/(3) _ p(3) + l(p(?) pW + p(l)p(2)s ) _ ^f^ 1 ) 3 

_(i?( 2 ) _ \F {1)2 ) U {1) - F (1) £/ (2) 

Again, the system of equations above can be consistently truncated anywhere. In fact the 
equation for only depends on fields of equal or lower perturbative order. Moreover the 
dependence on in the k— th equation is trivial. As already pointed out, i] only enters 
F^ 1 ', so that the first order explicitly depends on the noise, while higher orders are stochastic 
via the dependence on lower orders. Again, having adhered to an Euler scheme, results have 
to be extrapolated linearly as e — > 0. 

2.2.1. A numerical strategy. In the end NSPT simply amounts to integrate the equations 
(jll|) on a computer, i.e. we are now going to sketch the strategy for a perturbative Monte 
Carlo. If we want to perform a perturbative calculation to the order g n we need to replicate 
the gauge configuration n times, and evolve the whole set according to the Langevin system 
(jllj) . An important practical advantage of this method is that it can be coded by introducing 
very few modifications to a non perturbative Monte Carlo program for a local field theory. 
We have so far discussed the perturbative expansion of the Langevin equation. This choice 
is not essential: we could have chosen other stochastic dynamics, as long as they are based 
on differential equations (see for instance [T71 ITH1 ITT?]). Notice however that there is no 
possible perturbative expansion for an accept/reject Metropolis step, hence there is no NSPT 
analogous for such a kind of updating algorithm. In practice the mechanism of replicating the 
fields in perturbative orders consists in introducing in the field configuration the dependence 
on a new index (the perturbative order itself). Notice that any algebraic operation must be 
expanded perturbatively: 

X = A + B -> X {k) = A {k) + B {k \ 

k 

(12) X = A* B — > X (k) = aU) x B(k ~ j) - 

3=0 

Any measurement code for any observable 7 can be adapted to the present context by following 
the rules above. What we are typically interested in are the coefficients of the expansion 8 

(is) (o[J2g k 0(t)\)v = T,9 k o k (t), 

k k 

which also appear naturally as variables having a perturbative index and obeying the same 
multiplication rules (|12|). 

6 We move to a lighter notation and write for the Euler step U' = e~ F U. 
7 As for fermionic or gauge fixed observables see the relative sections. 
8 We adopt the notation for a generic field theory. 
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The average over the noise fields rj is reproduced, as usual, by averaging over a single long 
history, provided that one takes into account the autocorrelation at length n 

T 

lim(0 fc (*))„= lim l/TVOfcO'n). 

t— >oo T^oo ' 4 

i=i 

Finally notice that, if a non perturbative code is written by means of an object oriented 
language (where, typically, matrix operations are realized through the definition of suitable 
operators acting on classes), then the modifications just described could be really minor 
ones. 

A drawback of the method is the high request for memory (of course also the Floating 
Point operations requested by (fT2"j) increase with the perturbative order). These requirements 
however are achievable given our current computing power. 

The implementation of NSPT for fermions is quite different, and we will discuss it in 
section EJ 

2.2.2. The point of view of the algebra. The expansion (jl(Jj) may appear quite strange. As 
a matter of fact, one is quite familiar with a non perturbative formulation of Wilson action 
whose fields are the U XI1 , which take their values in the group. On the other side, Perturbation 
Theory amounts to taking into accounts fluctuations around a vacuum configuration, which 
naturally results in considering the Lie algebraic fields A XI1 . NSPT does not perform anything 
strange with this respect. One can also work with field variables which live in the algebra, 
where our perturbative expansion reads 9 

(14) ^(^-Er^M)- 

k=l 

This is perfectly equivalent. There is of course a perturbative relation between the two 
expansions: 

(15) A = tog(EO = log(l + X)0~* c/ ' ( * ) ) 

\ fc>o / 

3 

+ ... 

3 



The transformation (|15|) between the IPs and A's variables can of course be performed 
exactly at any finite order g n . We should stress once again the choice of the vacuum (back- 
ground) configuration which is embedded in (|15p. The Lie algebraic solution of the equations 
of motion is the usual (trivial) perturbative vacuum A Xfl = 0, so that the first non zero order 
is given by fluctuations of order 0(g) (i.e., in our preferred notation, 0(-^)). This, in turn, 

results in (fTUjl being expressed as 1 plus fluctuations of order 0(^=). It is worth to stress 
that in going from one notation to the other one has to handle the trivial series expansions 
of log(l + x) and exp(x), which can again be coded in an order by order notation according 
to the rules ()12j) . The constraint of unitarity on the original (not expanded) U xtl fields is 
translated into the usual (anti)hermitian, traceless nature of the A's fields 

A (k)\ = _ A (k) ^(k) = Q Wk 



9 Onc could observe that a better notation for the A field would be A xfl (t; 77) — ► J2k=i P~ k ^ 2 Ax/T (t] rj). 
This would enlighten the fact the first term in the expansion of A is actually a free field. Still, the bookkeeping 
of indices is more convenient in our notation, which is the reason for our choice. 
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i.e. the are still Lie algebraic fields. Notice that on the other hand there is no simple 
algebraic characterization of the fields, which are with this respect in a sense less 
fundamental. However, the transformation (jl5j) suffices to ensure that the fundamental 
relation UW = 1 = WU is satisfied at every finite order in the coupling constant (via 
highly non linear relations). Once expansions have been made, no matter what the choice 
of notation is, one has effectively decompactified the formulation of the theory (but this is 
not of course a feature of NSPT, but of Perturbation Theory itself). 

During the evolution it is necessary to periodically enforce the SU(3) group constraints, 
which may be spoiled by round-off errors. In our case this has to be done on the perturbative 
expansion of the gauge variables. The constraints are easier to enforce on the Lie algebra 
fields A^ k \ This is not the only advantage of working with variables living in the algebra. One 
can also get a non irrelevant saving of memory. In fact the highest order of the field typically 
only appears in observables linearly under trace, so that A^ kmax ^ can be omitted completely. 
On the other side the point of view of the algebra is less natural (in the end, Wilson action 
is formulated in terms of the f/'s). For instance the Langevin system of equations gets 
much more involved: 



(16) A®' = aM-FM 

A (2)' = A {2) _ p{2) _ 1 j-_p(l)^ 4(l)j 



2 

A (3)> = A (3)_ F (3)_1 [ F (l) jA (2)] _ 1 [ F (2)^(D] 

+^ [f« [F^UW]] +1 [A®, [fVa®]] 



For this reason we usually prefer the expansion entailed in Eq. even if for certain 

operations it is useful to switch to the algebra. 

3. Analysis of the stochastic process 

The discussion above shows that the perturbative coefficients in Eq. (fT3*|) correspond 
to well defined combinations of correlation functions of the stochastic processes 0^(t). We 
now want to stress once again that NSPT is a numerical implementation of Stochastic 
Perturbation Theory. With this respect not only the existence of limit distributions matters, 
but also the properties of convergence. In principle one would like to characterize as best 
as possible the stochastic processes in order to rely on NSPT as a computational tool. We 
will now restrict our attention to the case of gauge theories and in order to gain insight we 
will address in the following the study of the generic correlation function of M perturbative 
components of the fields 10 

M 

One can think of this expression as entering the perturbative expansion of a generic observ- 
able. Notice however that taken by itself the previous expression is by far more general. 
After showing that all such correlations have a finite limit for t — > oo, we will collect some 
results about the rate of convergence of these expressions. As a matter of fact, this will give 

10 In the following we will alternate the notations A xfi and A^(x). The latter is easier to read, in particular 
when there are other indices around (as in Xj) and it avoids confusion when Fourier transformations are 
involved. 
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us the chance to motivate some prescriptions in the implementation of NSPT for LGT, in 
particular the Stochastic Gauge Fixing, which is crucial, as already said, to demonstrate the 
results of The same questions were studied in ^T] for some simple models. The purpose 
of this section is to analyse the peculiarities of LGT. This will lead us to the introduction of 
Stochastic Gauge Fixing and the regularization of zero modes. 

Let us start considering the process defined by (fTTj) (or equivalently (fTBj)). with F given 
in (JBJ). It is not difficult to see that a limit distribution cannot exist. To illustrate this point 
we consider the system in a continuum regularization. The Langevin equation then reads 
(D^ b is the gauge covariant derivative; no expansion has yet been made) 

J^fax;*) = DfF^( V ,x;t)+r,;(x;t). 

The formal solution for the perturbative components of the fields (the solutions of the anal- 
ogous to (|TT|) and (jlEJ)) reads (in Fourier space) 

(18) A^ a (k; t) = T£ f ds e~ k2 ^f^\k, s) + I* f ds f^\k, s), 



p, \ "j fiu I Jv \ j I 1 p,v 

JO JO 

where Tj% and L ab v are the abelian transverse and longitudinal projectors 



k k 

rpab fx _ M v U rpabrpbc rpac 



k k 

jab C ^ X rab t be j ac rpab j be n 

fiv — _ p _0 «'" ^pv^vp — %p) I pu- L up — U ' 

The function f( n ' represents the interaction term, which only contains perturbative compo- 
nents of the field of order strictly lower than n 

f^ a (k;t) = gl^ n -^(k;t)+g 2 l^ n -^(k;t), 

f^ a (k;t) = r) v {k;t) a . 



Here Ip and 1^ are the n— th perturbative components of three and four gluons 
interaction. For instance 

gtf )a (k;t) = ^1 J dpd q 6(k + p + q )A b u (-p;t)AU-q;t)v^l(k,p, q ), 

where 

u uL<t(^)P> q) = fip.v{k — p) a + cyclic permutations. 
A first divergence appears projecting the formal solution ([18)1 by the longitudinal abelian 
projector L^ v . Along such directions the expanded Langevin system (j!6|) presents no damp- 
ing factor e~ fc2 * , for any of the perturbative components 

Lt^\k-t) = L% fdsf^\k,s). 

Jo 

The field A^ clearly diverges like a random walk (along such degrees of freedom). The 
behaviour of the higher perturbative components is difficult to predict in general. However, 
since the longitudinal components of A^ only depend on the {A^\m < n}, it is natural 
to expect diverging fluctuations at any order as was indeed observed in j^j. 

In the lattice regularization the interaction terms become much more complicated, but 
the argument for expecting a divergence remains unchanged. This problem is not avoided 
if we choose to work with the U variables (as in (jlljl ). instead of the algebra ones A (as in 
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([16))) . In fact, once the fields are perturbatively expanded, even the U variables do not live 
in a compact set anymore. 

A second source of divergence comes from the zero modes, since also for these no damping 
factor is present (see iJTBjl ). On a finite lattice the integration over momenta is substituted 
by a finite sum, and the degrees of freedom corresponding to k = give a finite contribution, 
enhanced in small lattices. This divergence results from a very similar mechanism as the one 
associated to the gauge degrees of freedom. 

None of these difficulties is peculiar of stochastic quantization. The first one appears 
whenever a perturbative expansion is introduced. In the traditional functional approach this 
implies the impossibility of defining the propagator, which is necessary to make perturbation 
theory meaningful. On a finite lattice, moreover, Perturbation Theory also faces the second 
problem: loops contributions are given by finite sums that in general entail singularities at 
k = 0. In both cases the phenomenon is related to the appearance of a zero eigenvalue in 
the action. In the stochastic approach a zero eigenvalue leaves the corresponding degrees 
of freedom without an attractive force, and diverging fluctuations (random walks like) can 
appear. This is not a problem if one is interested in the analytical computation of gauge 
invariant quantities. This was the spirit of the original proposal of Parisi and Wu [5]. 
However, what we have in mind is a numerical computation. From this point of view, 
diverging fluctuations - even in intermediate quantities - can be a serious problem, as it is 
clear from the plots in 

3.1. Stochastic gauge fixing. Stochastic gauge fixing was introduced in [T2j in order to 
provide the stochastic quantization approach with a non perturbative procedure of gauge 
fixing. The idea is to add a term to the Langevin equation such that the evolution of gauge 
invariant quantities is not affected. For instance in a continuum regularization we can write: 

= -SAj^- D ^ Vb[At]+ < {x '^ 

where F a [A,t] is an arbitrary non gauge invariant functional. The evolution of a generic 
functional F[A] is given by 

dF[A] f 5F[A] dA«(x;t) 
dt J 5A«(x;t) dt ' 

which is not affected by the presence of F a [v4,t] if F[A] is gauge invariant, that is 

One can also prove that the choice of Zwanziger 

(19) -Df V b = -Dfd u A b u 

* a ^ 

introduces a force that keeps limited the norms of the gauge fields. 

When we use a lattice regularization and the gauge variables live in the group, a generic 
gauge transformation has the form 

(20) U' xlt = e w * U XIM e- w *+A, 

where if is a field defined on the lattice sites and it takes values in the algebra. In this case 
a choice corresponding to ()19|) for the variation of the gauge field (modulo lattice artifacts) 
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(21) (0 < a < 1) Wx = - a J2 tfAvr 

The strategy in this case ^3] is to alternate a Langevin evolution step of order e with an 
order e gauge transformation chosen as (|21|). That means that a must be rescaled with e 
in the real simulations. Again one can show that such an operation introduces a force that 
keeps the norm of the gauge field limited. The iteration of (|2~T|) alone would drive the system 
towards a minimum of the gauge field norm 

N[w] = ^Tr (A^A^l) = ^Tr [(log U^) (log U^] . 

Notice that this results in fixing a popular gauge, since such a minimum is characterized by 
the Landau condition 

A 1 

There is a natural NSPT implementation of (J21j) : 

(o < a < i) Wx = £ /H™f w<F> = 

fc>0 m 

This prescription results in having also the gauge transformation implemented as a pertur- 
bative expansion. To be definite 

(22) e Wx = 1 + + . . . 

As expected, the (order by order) gauge transformation does not change the overall pertur- 
bative structure of the field, which is again 1 plus fluctuations of order 0((3~z). After the 
introduction of Stochastic Gauge Fixing the stochastic dynamics to implement is thus 

= e- F ^U x ,(n) 

(23) U Xfl (n + l) = e w ^U' Xfl e- w ^ u \ 

where one should keep in mind that everything is to be implemented order by order and F 
and w have to be taken as in Eq. (jHJ) and Eq 1)2 ip. Notice that there is a very simple way of 
looking at Eq (|2~3J) . which can be rewritten to 0(e) as U XfJL (n + 1) = e~ Fx ^ ua,Gr,G ^ U G (n), 
U x ^(n) being e Wx ^ U x ^(n) e~ Wx +^ u \ The effect of performing the gauge transformation 
after the Langevin step is equivalent to taking the Langevin step on a gauge equivalent 
configuration, but with a choice of rj given by e Wx rj x ^ e~ Wx (random field is gauge covariant, 
while equations of motion are gauge invariant). With this respect it is obvious why the gauge 
invariant quantities are not affected in the asymptotic limit of averaging over rj. 

Let us now proceed to understand what is the effect of Stochastic Gauge Fixing. This 
turns out to provide a force that keeps contained all the gauge degrees of freedom of all the 
perturbative components of the fields. In fact consider again how the Langevin equation (in 
Fourier space, in continuum regularization) is modified by the gauge fixing term: 

(24) A a (k;t) = -k 2 (T, u (k) + -L flu (k))A a u (k;t) + 

a 

+If a (k- 1) + t) + I^ a (k; t) + V ;(k- 1), 



Here d L is the backward derivative. 
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where ijf 1 (k; t) is the new three gluon interaction term introduced by gauge fixing. From 
the equation above one can obtain a system of equations as in and obtain a formal 
solution which reads (instead of (fTHj)) 



t ft 



(25) A< n \k-t) = T, u I dse- k2 ^f< n \k,s) + I ds e -£<*-) /«(»)(*, s). 

Jo Jo 

One can see that now the longitudinal (L^) degrees of freedom have a damping factor. We 
will use this observation in section I3~3~l where we will consider the problem of the convergence 
of the process. The reader is referred to the figures in [3] to have a feeling of how effectively 
Stochastic Gauge Fixing keeps fluctuations under control. 

The relation between Stochastic Gauge Fixing and the Faddeev-Popov prescription [20] 
was first studied in 21j. In section E~T1 we will show how it is possible to recover the case of 
the Landau gauge. 

3.2. Regularization of the zero modes. The problem of zero modes manifests itself in 
the formal solutions (fTSj) and (J2HJ- As already pointed out, the mechanism that leads to 
gradually diverging fluctuations is analogous to the one related to gauge degrees of freedom. 
Notice that in a finite lattice regularization the contribution of zero modes is finite, and 
particularly relevant for small lattices. It is interesting to notice that a similar problem 
occurs in finite volume lattice perturbation theory. The generic form for a given Feynman 
graph is in this case a finite sum which in general entails a singularity for zero momentum. 
A common prescription consists in simply dropping the k = contribution [22] ■ This 
is expected to reproduce the perturbative expansion of the lattice regularized functional 
integral in the infinite volume limit. In order to compute the perturbative expansion exactly 
associated to a lattice regularization in any finite volume one should define a consistent 
prescription in both perturbative and non perturbative context (for example with twisted 
boundary conditions 23j[23j). One possible choice for NSPT is to impose such a kind of 
boundary conditions. Another (simpler) approach is the subtraction of zero modes that can 
be obtained by imposing the condition: 

(26) JdxAf n \x) = 0. 

This means that there is no zero mode contribution to the field. This should be regarded as a 
prescription on a single mode, becoming irrelevant in the infinite volume limit for a quantity 
which does not have IR problems. This constraint must be imposed at each evolution step 
and for all perturbative orders. In fact the Gaussian noise rj^x) gives a contribution to 
the zero modes at each evolution step. Moreover the zero modes of higher orders take 
contributions also from non zero modes of the lower orders. For instance 

-t ig jabc 

'2(2tt) c 

n-l 



A« n \0;t) = I ds % ^-: + 







lb A n 

+ £/ dpdq5(0 + p + q)A b ^- 1 - m \-p ] t)A< m \-q ] t)v^ a (0, P ,q) + 

m=0 J 



+ 4gluons + G F. 
To keep condition ([26)1 . we enforce it after each evolution step. 



3.3. Convergence of the process. Now we come to the question of convergence of corre- 
lation functions of M perturbative components of the fields (|17j) . which in momentum space 
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read 

(27) ( n 4?>(M)>- 

1<3<M 

As already pointed out, the correlation function above is more general than the perturbative 
component of an observable, which is already known to converge ^2]. However these are 
the basic elements which build any observable, and we want to make sure that all of them 
converge, in order to have a safe numerical algorithm. Notice that we are not interested 
in the continuum or thermodynamic limits, at this stage. We are simply concerned with 
the existence of a limit distribution for every given finite lattice, where the simulations are 
performed. Therefore, for the rest of this section, a finite lattice is always understood. 

The proof of convergence was given in ^T] in the case of a 'zero dimensional' field theory. 
Since the argument is essentially the same, here we will only stress what is peculiar for a 
4-dimensional gauge field theory. 

First notice that any correlation function of free fields converges at least as 0(e~ q *), 
where q is the lowest Fourier component that contributes to the correlation. This result may 
be easily obtained by using the Langevin equation for free fields. 

Let us now define the total perturbative order of (j2*7|) as P to t = J2jPj- The idea is to 
reduce any function like ([27)1 to a sum of correlations of free fields. This is done by showing 
that (|27|) can be written as a sum of correlations that have either a strictly lower P tot , or the 
same P to t but a strictly lower number of fields. 

First of all we write (|25|) in a discretized Langevin time t = Ne (and a = 1). 

(28) Af{k-t) = e- k2 *A^(k;t-e) + ^%(k;t), 

Af(k-t) = e- k ^A^(k;t-e) + ef^(k;t). 



It is worth noting that the previous expression is valid for all degrees of freedom (and with 
all k 2 > 0) only because of the gauge fixing and zero momenta subtraction described in the 
previous sections. 

By inserting (|2*Hj) in (}2Tj) and taking first the limit e — > (at t — Ne fixed) and then the 
limit t — > oo one finds: 



(n^i)> 



(29) 



Z^m=l h m 



E E < II < J (W) + 



{l<i<M\pi=0} {i<j<M\pj=0} {l<h<M\h^i,j} 



e < n ^\k h )f^\k h )) 

.{l<j<M\ Pj ^0} {l<h<M \h^j} 



The argument to deduce (|29J) is given in [TT| . for the case of a 0-dimensional field theory. 
The only difference here is that fields are now labeled by discrete momenta (remember that 
we are in a finite lattice). As a consequence the correlation functions of two fields have a 
decay constant which depends on momenta and gives the momentum dependent factor in 
front of the square bracket. 

By iteration of (|29*j) Ptot times, the initial correlation function is reduced to a finite sum 
of correlation functions of free fields. This proves that any correlation function ()27jl has a 
finite limit for t — » oo. 
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We now come to the point of the convergence rate. For free fields the convergence to 
equilibrium may be checked directly and it turns out to be 



(JJ A^( Xj ; t)) = floo + a ie - q2t + a J e ~ 



where q 2 is the lowest Fourier component entering the correlation function. Correlations in- 
volving higher perturbative components of the fields also have the same exponential damping 
factor. However also power corrections appear, such as 

t p e~ q2t , 

where p is of the order of P tot . 

For a given observable at a given perturbative order, it is therefore possible to estimate 
the time scale in which convergence to equilibrium should occur. This is in practice very 
important. Of course one would like to know the size of all the moments of the distributions, 
which tells the size of the fluctuations at equilibrium. However the determination of such 
moments is not easier than determining the perturbative coefficients themselves, and even 
for the toy models studied in [TT] only partial results could be attained. The analysis of 
the fluctuations must be done for any new observable, by means of the statistical analysis 
described in HTJ. 



We want to stress the reassuring message coming from llj: one does not observe in 
the case of LGT the pathological behaviour of very large (and not normally distributed) 
fluctuations at high orders which we found in very simple (0-dim) models. In the spirit of 
illustrating what are the typical fluctuations one has to live with, in the following we perform 
an error analysis on a benchmark computation, i. e. a large set of Wilson loops. 

3.4. Analysis of performances and simulation time. Before analyzing the fluctuations 
of a particular observable, let us consider the computer time needed for a single Langevin 
iteration. Together with autocorrelation this will provide a precise estimate of the cost of 
reaching a given precision. A quick inspection of the algorithm suggests a dependence on 
the linear lattice size (L) and the maximum perturbative order (p) as 

(30) isKi4 .(P!^l. 

In fact the algorithm for pure gauge fields scales with the volume of the lattice, and most 
of the time is spent in performing order by order multiplications. This formula fits well the 
timings we measured (see Table Q). The above considerations are of course not conclusive, 
since they only deal with crude execution times. In order to further comment on the efficiency 
of the method, we present in the following the error analysis of a large set of Wilson loops. 
In this context we will also add something about autocorrelation. 

3.5. Wilson loops in quenched NSPT. The results on Wilson loops we are going to 
present were used in [§] to compute in Lattice Perturbation Theory the static potential, 
from which the static self-energy was in turn computed. The results for the Wilson loops 
themselves have not yet been published till now. Apart from the interest for reconstructing 
the static potential, Wilson loops have many other reasons to be regarded as interesting (for 
a discussion and a good list of references see for example [22] )• On top of that they are in 
a sense also a typical benchmark. One can find other perturbative computations of these 
quantities, either in the standard approach (even an incomplete list of reference should at 
least cite [2511221123122]) and by a less conventional technique Q28J). To our knowledge the 
list of results we present here is the largest available at the moment. 
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The expansion (jlUj) was made up to the maximal order k = 6, so that we could obtain results 
up to 0(/3~ 3 ). This was done on a 32 4 lattice. This size was the largest which was possible to 
simulate at this given order (third, as counted in the physical expansion parameter) on the 
computing facility available to us at that time (it was an APE100 in the crate configuration, 
with 128 Floating Point Units (FPU's) and extended memory, for a total peak performance 
of 6.4 x 10 9 Floating Point operations per second, i.e. 6.4 GFlops). Notice however that any 
execution time we report in this paper refers to the APEmille architecture, which is the one 
currently in use. 

Measurements are taken by following a common procedure in Monte Carlo simulations. 
One first lets the system thermalize and then starts collecting configurations which are stored 
and used to measure various observables. Notice that in the case of NSPT this results in a 
database which asks for a big storage area. Since the fields are replicated in perturbative 
orders, one single configuration of a 32 4 lattice up to sixth order amounts to 1.7 Gbytes of 
binary data 12 . As already pointed out several times, configurations have to be collected at 
least at two different values of the Euler time step. An important issue is of course the deci- 
sion on the frequency at which one should write out the configurations. We usually measure 
Wilson loops (for which many results are already known at least up to g A order) to decide 
when the system is effectively thermalized (one can check several sizes and so several effective 
scales). We also monitor on the fly at least the basic plaquette (which is anyway measured 
to compute the equations of motion entailed in F) in order to estimate the autocorrelation 
time at least for such a (very) local quantity. Some autocorrelation measurements for the 
plaquette on a 32 4 lattice are given in Table El We did not attempt a precise determination 
of the autocorrelation of each Wilson loop, since no long enough history is available. Instead 
we analysed them by the bootstrap method. NSPT results should in general be analysed 
by means of the bootstrap method as described in jTTj. This is necessary also to take into 
account the possibility that fluctuations are not normally distributed. As already said, how- 
ever, by inspection of frequencies histograms up to a 10 it turns out that no huge deviations 
from normality are observed in the case of LGT. In fact the bootstrap analysis coincides 
with the traditional one. This phenomenon is not difficult to understand if one considers the 
origin of the large fluctuations in O-dimensional models [TTj . 

We collect our results for all the Wilson loops up to 16 x 16 in Appendix A. Of course 
there are larger and larger finite sizes effects. Comparing results with smaller lattice sizes 
suggests that for loops up to the extension 6x6 finite size effects are roughly of the order 
of the statistical errors. Deviations of the order of a few percent occur at g e for both loop 
sizes of the order of half the lattice size. Results were obtained out of measurements of 120 
configurations. By inspection of the results one can see that relative errors range from order 
10~ 4 for leading order up to 1CT 3 for third order. Actually the ratios of relative errors at 
different orders stay much the same for various loop sizes. At fixed order, relative errors 
increase roughly linearly with the loop size (perimeter, actually), which is quite reasonable. 

4. Other applications 

4.1. Perturbative expansion for gauge-fixed observables. Once Stochastic Gauge Fix- 
ing has entered our prescriptions for NSPT simulations, any correlation of gauge fields con- 
verges and also non gauge invariant quantities are calculable. It is interesting to inspect 
how one gets a stable signal also for non gauge invariant quantities: in figure Q we give 



; The space requested can be substantially reduced (a factor 8/18) if one stores the Lie algebraic fields. 
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Figure 1. History of the gauge dependent quantity TrU^, order (3 (4 lattice). 



an example of the evolution of a gauge dependent observable (the trace of the link). Even 
if it is not easy in general to make contact with the standard covariant gauges result, the 
latter can be easily recovered in the case of Landau gauge. The strategy is the same used 
in standard (non perturbative) Monte Carlo simulations. One first produces a thermalized 
configuration and then fixes the gauge. The Landau gauge condition is easy to attain if 
one consider Eq. (|2ip. If iterated, that gauge transformation drives the minimization of the 
norm functional, which results in turn in the Landau condition 13 . This gauge transformation 
is by far more effective if Fourier accelerated, as pointed out in [THj. Figure El shows how 

the quantities y J2 X Tr {(J^^^fl^^Y (X^^u^m)) are minimized by the iteration of (the 
order by order version of) Eq. (|21jl . ensuring to machine precision the Landau condition at 
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Notice that this is true at every perturbative order. 
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every perturbative order. 

We make a last point on gauge fixed computations. As it was first proposed in a very 
appealing feature of NSPT is that there is also a way to obtain covariant gauges results. 
Without entering into details, it suffices to recall the basic formula of the Faddeev-Popov 
procedure [20] , i.e. the complete measure in the path integral 

e -(s G +s gf ) detL[E7] 

where one usually chooses (covariant gauges) S g f = ^ X^s^^A^) 2 - m the previous for- 
mula a gauge fixing action has been added to the gauge (in our case, Wilson) action, but this, 
as we know, is not the end of the story: one must also keep into account the determinant 
of the Faddeev-Popov operator. Even if one can always define a new contribution to the 
action by 

detL = e- ( - TrlnL) 

the previous expression is in standard Perturbation Theory almost useless. One trades the 
non locality of this action with the introduction of the (spurious) ghost degrees of freedom. 
In NSPT one instead makes direct use of the (Tr InL) action. Since the mechanism is just 
the same that makes it possible to treat also the fermionic determinant, we differ a few 
comments on this point to section 5. 

4.2. Perturbative expansion around non trivial vacua. Any perturbative expansion 
must be performed around a fixed vacuum configuration. One nice feature of NSPT is 
its flexibility with respect to different choices of the perturbative vacuum. In this very 
brief section we will show how, from a technical point of view, NSPT only needs very few 
modifications when the perturbative vacuum is not the usual one (defined by U = 1). This 
is not difficult to understand, since a vacuum configuration is basically a solution of the 
equations of motion around which the solution of Langevin equation fluctuates in force of 
the random noise. In order to illustrate this possibility we will just give the relevant recipe 
in one example: the vacuum configuration of the Schroedinger Functional scheme [HO] . 

In the approach of (3U1 the following boundary conditions are fixed for links pointing in 
all space-like directions at the (boundary) time slices (x° = and x° = T) 

U(x,o)k = exp(C), £%,T)fc = exp(C"), k = 1, 2, 3 Vf . 

We are not interested here in the precise form of the C and C (which are £77(3) matrices). 
The classical solution with these boundary conditions is 

(31) V x0 = 1; V xk = exp([x°C" + (T - x°)C]). 

The strategy in this case simply consists in replacing (fTU|) with: 

U Xfl (t; V ) - V x , + 0- k ' 2 UW(t;ri). 

k=l 

Apart from that, only slight modifications are needed. We are not going into details here. 
Figure 01 shows a typical signal for Schroedinger Functional NSPT (it is the leading order 
plaquette). We conclude with a last, trivial remark: it makes no difference how complicated 
the vacuum configuration is. Such configuration could even be known only numerically, for 
example from a quenching procedure minimizing the equations of motion. 
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Figure 3. First non trivial order (J3 x ) of the plaquette in presence of the 
vacuum configuration defined in (|3T|) (8 4 lattice). 



5. Inclusion of dynamical fermions 

We now come to the discussion of how to include the contribution coming from fermionic 
loops. This was first discussed in [31 and a first, preliminary discussion of the implemen- 
tation can be found in [32J. Let us first of all define our notation. We write M[U] for the 
fermion matrix (i.e. the Dirac operator). Even if most of what follows is valid for a wide 
class of fermionic action, we will specify to the case of Wilson fermions (which is the case 
for which we will present prototype results later) with Wilson parameter r = 1: 



Xfl 



Let us also write down the explicit form of the matrix element of M[Z7]: 
(32) M yf3biZJC [U] = ^2(m + 4)5 ys 5 fh 5 bc + 

X 

4 

-h/, {$V,z+i> i 1 + lv)Pl ( U L)bc + 5y,z-D (1 - 7v)/3 7 {U Z -uu)bc) ■ 
v=\ 

In the following we will drop in our notation the dependence of M on U. We will also make 
no comment on the dependence on the number of flavors rif, since it is trivial to work it 
out. Actually most of the NSPT simulations that we have been running involve a couple of 
degenerate massless quarks. This is the most relevant case, for example for the computation 
of renormalization constants for QCD. The functional integral for the fermionic degrees of 
freedom can be computed, resulting in the well known fermionic determinant contribution 
det M. It is worth to recall here the point that we have already made in the previous section: 
the Faddeev — Popov action has the same determinant structure and this means that 
the strategy for including fermions is just the same needed to implement the Faddeev-Popov 
procedure (without ghost fields). We will come back to this point later. 
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Once the fermionic contribution has been integrated out, the result appears as a determi- 
nant and one has to manage a path integral which is integrated only on the gauge degrees 
of freedom, but with a new weight given by 

(33) e- SG det M = e~ s ^ = e -(Sa-TrlnM) 

In the previous formula one has rewritten the determinant as a contribution to a new (ef- 
fective) gluonic action. Let us see what is the effect of this on the prescription for writing 
down the Langevin equation (0). One simply needs to replace 

V"Sg^V»A// = V*S G -V*TrlnM 



xm°G m V^Oe// - V xfl O G - V aAe 

7 x^G — Tr (( . r/r 



(34) = V- Se-Tr ((V" M)M~ l ) 



Eq. ()34p asks for computing the trace of a product of two matrices: the Lie derivative of 
the Dirac operator and the inverse of the Dirac operator itself. The real difficulty does not 
come from the Lie derivative, which is pretty simple, since it is substantially local: 



i 

2 

&y,z—ji &x,z— fi (1 T/n)/37 (-^ ^Xfjbjbc , 



(35) VyiW 7C = % - ( 8 yjZ+ii 8 z , x (l + ltl )p 7 (Ul lx T a ) bc + 



The difficult task comes from having to face an inverse, i.e. non-locality. Eq. (|54j) and 
the solution given to it by the Cornell group fH] were in a sense a prototype for fermionic 
simulations. The idea is to introduce another Gaussian source £ 

which enters the new version of Eq. (JHJ) , which reads 

F = T a (e$ a + yfeq a ) 

(36) <£ a = [V a xp S G - Re (^(V^MMM-^Un)} . 

In the previous formula all the repeated indices are to be summed over, keeping also in mind 
that k,l,n are multi-indices, to be intended as in Eq. (|32j) . This of course means that the 
field £ has got, on top of position (or momentum), the degrees of freedom of a so-called 
spin-color field. The evolution of the process will now average both on rj and on £. The 
average on £ will in particular results in 

= [VlpSc - Tr ((V^M)M- 1 )] 
= V^[S G -Tr(lnM)] 

in which one can recognize Eq. (|34|) . One can now rewrite Eq. (|3*Bj) as 

(37) $ a = [Vl^Sc - Re (^(V^M)^)] 
where the vector if> is the solution of a linear system 

As already pointed out, this was in a sense a prototype solution for fermionic simulations of 
Wilson fermions action. One has ended up with the inversion of a sparse matrix on a given 
vector. Of course the fact that the matrix M is sparse is crucial for the actual efficiency of 
the method. 

We now proceed to our NSPT version of the algorithm. Both Eq. (J34j) and the prescription 
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entailed in Eq. (j3T?j) can be simply translated according to our order by order prescription. 
This means first of all that the matrix M itself gets expanded as a power series 

(38) M = M (0) + P~ k/2 M {k) . 

k>0 

By direct inspection of Eq. ()32|) one realizes that the expansion of M is trivial (it depends 
linearly on the U xtl ): 

x 

4 

-^^(^(l + ^^^t). + S y ^(l- lu )^(U { z % u ) bc ) . 

The only non-trivial feature of the previous expression is that also the mass acquires an 
expansion. This is due to counterterms, as we will discuss later. Much the same holds for 
V^M, as it appears from Eq. (}35|) . One now has to face the inversion of M as a matrix 
power expansion. It is easy to compute 

m- 1 = ^/r fc/2 M- l(fc) 

(39) = M^ 1 + J2^ k/2 M~ l(k) . 

k>0 

The notation enlightens the fact that the zeroth-order of the inverse is the inverse of the 
zeroth-order. As for higher orders, a simple recursive relation holds: 

M -i« = _ M (o) 

M -l(2) = _ M (0) 

M" l(3) = -M® 

(40) M" l(n) = -M(°) 



We now mimic the construction in Eq. (j3*Uj) introducing the new random field £. Notice that 
this is a field with no power expansion. Once the product (V"„M) (M -1 ) has got translated 
in a matrix power expansion (i.e. it has become a matrix sum of perturbative orders), the 
£'s are in charge of taking the (stochastically evaluated) trace of the various orders, which 
results in getting the power expansion of Tr ((V" At M)M~ 1 ). Just as in Eq. (J37)l it was useful 
to introduce the vector ip, it is now useful to define a (power expanded) new field 



- l M^M- l{1) 



n-l 



1 M (n ~ j) M U) 

j=0 



-1 



(41) 
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We now have a compact expression for the n th order of the fermionic contribution in Eq. (JHTj) 14 . 
which reads 

n 
3=0 

Having already made the point that a simple recursive formula holds for Eq. (J4Uj) . it is 
straightforward to notice that a similar relation holds for the computation of the fields ift^); 



(42) 



^ 2 



if, 



(n 



M (o)-^ 



ra-1 



3=0 



5.1. The actual implementation of unquenched NSPT. The structure of the original 
(non perturbative) Eq. (J37j) is that of a scalar product: one starts constructing ip (which is 
the solution of a sparse system), then apply the matrix (V" M) to it and finally contracts 
with £. The NSPT version has of course the same overall structure, provided that every 
operation is intended as order by order (and with the caveat that the field £ has only zero 
order). We now proceed to illustrate why there is a natural and efficient way of implementing 
this program on a computer. 



• Eq. (|42|) (or Eq. (|40)l ) makes it clear that there is no matrix inversion to undertake. 
The only inverse matrix is the zero order 1 , which does not depend on the fields 
and is the standard tree level Feynman propagator, diagonal in momentum and color 
space 

_ (m+fk 2 ) -i)jt 
(m+fk 2 ) 2 + k 2 

The only problem associated to this expression is its computation at k = in case 
of zero bare mass. Various infrared regularizations are available: leave a small mass, 
introduce anti-periodic boundary conditions for the fermionic fields in time direction 
or constrain the zero mode degree of freedom to zero. To produce the results in the 
present work the last choice was adopted. We checked stability of results with respect 
to the use of a small mass. More on this in section 6.2 

• Notice that Eq. (J4*2~j) suggests the construction of the various orders ip^ in a sequen- 
tial way. At every order only one application of is needed and the propagator 
operates on a sum of already computed quantities (i.e. lower order i/j's). While 

M(°) is diagonal in momentum space, all the operators needed to construct this 
sum (i.e. the various orders of M) are almost diagonal in configuration space. This 
suggest the strategy of going back and forth from Fourier space. 



k and I are again dummy (multi-)indices over which a summation is understood. 
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• This structure is common to many fermionic regularizations. This is also the right 
point to go back to the Faddeev-Popov determinant. In the latter case the zero order 
inverse is again reminiscent of a tree level propagator, but this time a scalar one 15 . 
As for higher orders, these are in the case of Faddeev-Popov determinant naturally 
expressed in the adjoint representation. The strategy of going back and forth from 
Fourier space stays much the same. 

• Previous considerations clearly stress the need for a reasonably efficient Fast Fourier 
Transform (FFT). Apart from the stochastic dynamics, actually momentum space is 
the natural stage also for the measurement of observables which are diagonal in that 
space (think about quark bilinears). 

• By inspecting Eq. (|36|) one can not recognize whether the effective time scales asso- 
ciated to the gauge and (pseudo)fermions drifts are the same. This is a very general 
and well known point [33]. Besides this observation, even if one decides to make use 
of the same time step for both drifts, it is anyway useful to take advantage of the fact 
that one can always live with 0(e) errors. In our Euler scheme this is an effect which 
has to be in any case extrapolated to zero. As a matter of fact the computations of 
the gauge and (pseudo)fermions contribution to the drift are quite different from the 
point of view of the actual implementation on a computer. This suggests to break 
down the evolution step in the following way: 

— Evolution by the pure gauge contribution to the drift F gauge = e VSg + \/^V- 

— Construction of the ip^; this is the only non-local piece of the computation, 
the non-locality being anyway traded for multiple applications of an FFT, after 
which every operation is indeed local. 

— Evolution by the fermionic contribution to the drift Ff erm = eRe(^(VM)i/)). 
This does not present any structural difference with respect to the first module. 

One should always keep in mind that also a step of Stochastic Gauge Fixing has to 
be taken at the end of this sequence. 

There is a last point to be made concerning the NSPT expansion of the unquenched Langevin 
equation. In inspecting the structure of 

F = e(VS G - VTr InM) + Jli] 

one should keep in mind our rescaling of the time step e' = e/3 which was meant to get a 
consistent perturbative expansion. Notice that this leaves the fermionic contribution to the 
drift with an overall f3~ l in front. This does not come as a surprise. Once the fermionic 
degrees of freedom have been integrated, one is left only with the gauge bosons. In the 
equation of motion for the gluons fermionic contributions enter as loops and to "dress" a 
gluon one needs at least an order 

5.2. Analysis of simulation times. We saw in Sec. 3.4 that the simulation time for the 
pure gauge implementation of NSPT is in good agreement with the expectations: it scales 



This does not come as a surprise: the determinant involves the degrees of freedom which in standard 
perturbative computations are absorbed by ghosts. 



23 



Table 1. Simulation times in (seconds/iteration x number of processors). 
Run in 8 4 lattice were performed on an APEmille board; 16 4 on an APEmille 
unit; 32 4 on an APEmille crate. Absolute values in seconds are only for il- 
lustration, being strongly implementation dependent. One should observe in 
particular the scaling in L, in p and compare Quenched with Unquenched. 



lattice size L 


order p 


n f = 




8 


9* 


27.8 


48.0 


8 


9 8 


49.0 


81.6 


8 


9 W 


78.7 


128.7 


16 


9" 


453 


814 


32 


9" 


7168 


12979 



with the volume and is dominated by the order by order multiplications (see Eq. (j3*Uj) ). 
Since the unquenched version of the method introduces substantial changes, it is of course 
compelling to verify what is the impact on the simulation times. 

This is a good point to comment on our own implementation of NSPT programs. The main 
NSPT project for Lattice QCD is run on the APE architecture. The first implementation 
was that of the quenched version on the APE 100 family. The unquenched version has been 
developed on APEmille. Our FFT implementation mimic |33] . which is based on a 1-dim 
FFT plus transpositions. The latter operation is the one asking for local addressing on a 
parallel architecture, which made it necessary to wait for APEmille in order to implement 
unquenched NSPT on APE. collects some other comments on our APEmille implemen- 
tation. In the last two years the growth in the computational power available on PC's made 
it worth to develop a C ++ implementation which is now run on medium size PC-clusters, 
usually to assess finite volume effects (we simulate small lattices on PC's and large lattices 
on APE). In Tabled we make use of timings taken on ^4PP in order to assess what is the 
overhead in moving from quenched to unquenched NSPT 16 . We report the execution times 
for a single iteration times the number of processors (formally, a theoretical execution time 
on a single processor). We are mainly concerned in scaling properties. On a given volume 
(which in the example of the table is a modest 8 4 ), one can inspect the growth of execution 
times as the perturbative order grows. Both in the column of the quenched and in the col- 
umn of the unquenched simulations, the scaling of computational time is again consistent 
with the fact that order by order multiplications are the dominant operations. On each row 
this results in a growth in time due to unquenching which is roughly consistent with a factor 
5/3. One then wants to understand the dependence on the volume, which is the critical one, 
given the presence of the determinant: this is exactly the growth which has to be tamed 
by the FFT. One then compares execution times at a given order on L = 8, L = 16 and 
L = 32 lattice sizes. Notice the different APEmille configurations: L = 8 is simulated on an 
APEmille board (for a total of 8 FPU's), while L = 16 on an APEmille unit (32 FPU's) and 
L = 32 on an APEmille crate (128 FPU's). One easily understands that FFT is doing its 
job: the simulation time scales with the volume also for the unquenched version of NSPT. 
The very same emerges also from the (C ++ ) PC-cluster implementation of the programs. 

As already pointed out for the quenched case, at this level one has only compared crude 
execution times. In our experience the optimization of the signal to noise ratio asks for 
smaller values of the Euler time step. We halve the values of time step with respect to the 



'The quarks involved were degenerate in mass, so that the dependence on n/ is trivial. 
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Table 2. Autocorrelation times for the basic plaquette expressed in iteration 
number. 



lattice size L 


n f 


Euler time step 


r 1 


P~ 2 


r 3 


32 





0.01 


~ 50 


~ 70 


~ 100 


32 


2 


0.005 


~ 100 


~ 120 


~ 150 



quenched case, which for the computational overhead results in an overall factor 2 on top 
of the 5/3 coming from the crude execution times measurements. This is anyway a good 
message: we are talking of a difference with respect to the quenched case which is a factor 
and not an order of magnitude. Table 121 contains estimates of autocorrelation times for the 
basic plaquette both in quenched and unquenched case. If one rescales the values keeping 
into account the time step, the unquenched autocorrelation time appears slightly shorter. 
Even if this is not conclusive, it is not unreasonable, since the new random field £ can give 
an extra contribution to decorrelate two successive configurations. 

A more complete assessment will be possible after having inspected the errors in a bench- 
mark computation which is the same also reported for the quenched case, i.e. the computa- 
tion of Wilson loops. 

6. Results in unquenched NSPT 

We collect, in the following, a couple of prototype applications of unquenched NSPT: the 
computations of unquenched Wilson loops and of the (Wilson fermions) critical mass (to 
be introduced later). Both are presented to order /3~ 3 for the case of Wilson gauge action 
coupled to two mass degenerate Wilson fermions. The quark masses are put to zero by 
plugging in the relevant counterterms (more on this in the following). Actual computations 
were performed on a 32 4 lattice on an APEmille machine in the crate configuration (128 
processors for a peak performance of 64 GFlops). 

6.1. Unquenched (jif = 2) Wilson loops. In Appendix B we report a table of Wilson 
loops which are just the same as in Appendix A, but in the unquenched case of two massless 
Wilson quarks. As in the quenched case, the major physical motivation was the compu- 
tation of the static potential and the static self-energy [22] • The configurations on which 
measurements were taken were 200 17 (again, at two different values for the Euler time step). 
By inspection, the general picture for relative errors does not change with respect to the 
quenched case. 

6.2. The ri/ = 2 critical mass for Wilson fermions to the third loop. It is well known 
that the Wilson lattice regularization of fermions breaks chiral symmetry. This comes from 
the irrelevant term which enters the Dirac operator in order to cure the so called doubling 
problem. The first net effect is then the appearance of an additive mass renormalization, 
which is usually referred to as the critical mass. 

Let us set up our notations. We write the two points vertex function (the inverse of the 
quark propagator) as 

r 2 (p 2 ,m) = S{p 2 ,m)~ l 
(43) = i $ + m — S(p 2 , m) 



17 These configurations are part of a wider database containing also other values of n/. 
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where 

(44) S(p 2 , m) = S c + m Y> s (p 2 , m) + i^T lV (p 2 , m). 

The previous equations are written in the continuum limit. S c is our notation for the crit- 
ical mass. In our simulations m = 0, so that S c is the only contribution along the (Dirac) 
identity operator. For a ^ 0, one has 18 S c i— ► S c + S c (pa). By restoring physical dimensions 
one can inspect the a~ l divergence in S c . Because of the power divergence, the perturbative 
evaluation of the critical mass is not supposed to be good. Therefore this quantity was in 
a sense a prototype for non perturbative determination of renormalization constants (an 
additive renormalization, in this case). Still, its perturbative computation has become sort 
of a benchmark. Two loops results are available [HBJ EZj so that we have the chance both to 
check the first and second order and to try to understand what is the effect of the third loop 
on the (poor) convergence properties of the series. 

The computation is performed on sets ranging from 200 to 60 configurations, depending 
on the value of momentum 19 . These are part of the same configurations database on which 
also Wilson loops were measured. As already said, the mass of the two quarks was put to 
zero by plugging in the relevant counterterms. This is a good point to comment on what this 
means. Once one realizes that there is an additive renormalization entering the stage at one 
loop, each new result for this additive contribution to the mass has to be taken into account 
in the computation of higher loops. In conventional, diagrammatic Perturbation Theory this 
means to take into account new effective vertices (counterterms, actually). NSPT is not dif- 
ferent with this respect. The higher the loop order we want to go, the more counterterms we 
have to plug in coming from lower loops. Since we know the first and second counterterms 
from [SniEZl, there is no problem in trying to compute the third loop. Also the latter should 
in turn be plugged in if one wanted to compute the fourth order correction. Notice that 
plugging the counterterms in results in what one would call a renormalized perturbation 
theory (with respect to the additive mass renormalization). Since counterterms are taken 
into account, the result that we obtain for the first and second order of the critical mass is 
zero, as can be seen from figure |U This is the way to check that we agree with previous 
computations at first and second loop level. 

In order to compute S c we first compute the propagator and then invert it to obtain IV 
From T2 the critical mass is obtained by taking the projection along the (Dirac) identity 
operator. The computation is performed in the Landau gauge by using the gauge fixing pro- 
cedure that was explained in section 14.11 The computation of a propagator is the prototype 
measurement in the fermionic sector of QCD. We know from the Wick theorem that every 
observable reduces to a combination of inverses of the Dirac operator. This is the case in 
which the measurement we want is the propagator itself. In NSPT we measure it in mo- 
mentum space. Since the configuration average will be diagonal in momentum, it suffices to 
compute on each configuration only diagonal entries of the inverse of M. Of course we still 
make use of our recursive relations given in Eq. ()40j) to compute the inverse of M. Actually, 

18 Also the tree level mass for a ^ has an irrelevant contribution due to the Wilson prescription to 
eliminate doublers: m i— > m + mw(ap). To shorten the notation in the following we will write £ c (pa) = 
S c + £ c (pa), i.e. S c will be both the function of pa and its value at pa = 0, which is the relevant, although 
divergent, physical quantity. 

19 This means that we have not evaluated the propagator on all the 200 configurations for all the momenta 
that are plotted in figures 0] and [5] As it can be inspected, this does not appear manifest in the end, because 
the quality of errors is quite good. 
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Figure 4. First and second loop of S c after the relevant counterterms have 
been plugged in. Diamonds are the fitted points. Error bars on the measures 
are hard to distinguish. 



we make use of the recursion given in Eq. ([42)1 . since the various entries of the inverse matrix 
are obtained by operating on sources vectors £ which are now (Kroneker) (^-functions. This 
is the trivial observation that the Ay entry of the matrix A can be computed as a scalar 
product 

k,l 

where the entries of the vector are defined as ^ = 5 a k- We make use of the symmetry 
properties of the propagator by averaging over the momenta that are connected by hyper- 
cubic group transformations. This is not the only reason why the hypercubic symmetry is 
important. As a matter of fact, in the continuum limit one recovers the 0(4) Euclidean 
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Figure 5. The third order critical mass T,^\ As in Figure EJ diamonds are 
the fitted points (error bars are again very tiny). The point marked as a cross 
is the extrapolated value. 
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symmetry, i.e. £y an d 51 s are functions of p 2 (logarithms, actually). At finite values of the 
lattice spacing only hypercubic symmetry is there and as a consequence the p-expansion 
for a generic entry of the propagator contains components on all the hypercubic invariants 
one can construct from powers of p. By restoring physical dimensions one realizes that the 
expansions are in power of ap. This means that we have a handle on lattice spacing effects 
and this is the way to look at figures 0] and 03 We write the expansion of the critical mass 
as 

(45) -e c = E^/r 1 + ^ 2) /r 2 + sf/r 3 + • • • • 

Once again, what we plot in the figures are not the coefficients of this expansion, but those of 
— S c — £c 1 ' ) /3~ 1 — [3~ 2 . These coefficients are zero at first and second order (see figure HJ), 
while Sc 3 ^ is the first non zero coefficient, i.e. the one we are interested in (see figure EJ). No- 
tice that these coefficients are plotted versus (ap) 2 , but, due to the presence of higher orders 
invariants, they do not appear smooth. Since we are interested in the continuum limit of the 
computation, what we have to do is to extrapolate the curves to 0. This is done by fitting 
the data taking into account the hypercubic invariants made of powers of ap. The curves 
in figure 0] actually extrapolates to numbers of order 10~ 3 and 10~ 2 . We explicitly checked 
at one loop level that within errors we exactly reproduce the finite size effects expected on 
a 32 4 lattice. The goodness of the extrapolation to zero of the first two coefficients (with 
counterterms taken from an infinite volume computation [3UJIIIZ]) is indeed a good indication 
on finite volume effects and sensitivity to our choice of infrared regularization (see section 5). 

The final result for the third loop of the critical mass for n/ = 2 is 

(46) = 11.79^;. 

The errors are estimated by looking at the stability of the fitted result with respect to varying 
the number of points included in the fit and the higher orders lattice invariants taken into 
account. They are only statistical errors, but the finite size effects are expected to be small. 

In order to inspect the convergence properties of the series (which, as we have already said, 
are supposed to be poor) one could now try to repeat the analysis contained in J361 , ^. e. the 
comparison of a fixed order in Perturbation Theory with the non perturbative determination 
of the critical mass. If one for examples tries to perform this little exercise at (3 = 5.6 using 
for the non perturbative determination of the critical mass the data in [SE], the ratio of 
the perturbative result to the non perturbative result raises from 0.55 at first loop to 0.71 
at second loop and finally to 0.79 at third loop level. As one knows, results can improve 
quite a lot moving to boosted perturbation theory. We do not perform this exercise for 
the critical mass and leave it for the reader, wanting to stress that this is not physically 
relevant in this case. Still, we want to stress that what we have just presented is only a 
prototype computation. Actually, by considering the projection of T 2 along the 7-matrices, 
one has access to the field renormalization constant. Preliminary results have been reported 
in (HHI and the complete computation of the latter will be reported together with that of 
renormalization constants for quark bilinears. For those logarithmically divergent quantities 
a three loop computation would be certainly useful. 

7. Conclusions 

We reviewed the application of Numerical Stochastic Perturbation Theory to Lattice 
Gauge Theories and discussed the implementation of the method in a Monte Carlo algo- 
rithm. We studied the underlying stochastic process and its properties of convergence. In 
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order to assess the efficiency of the method we measured crude execution times and studied 
the size of the fluctuations for benchmark computations. After touching the subject of gauge 
fixed computations and expansions around non trivial vacua, we discussed the inclusion of 
dynamical fermions. We stressed that in this perturbative approach they are not so computa- 
tionally expensive as in the non-perturbative case. We presented as benchmark unquenched 
computations a large set of Wilson loops and the critical mass of Wilson fermions to the 
third loop. 
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Appendix A. Quenched Wilson Loops 



Loop dimensions 


order (3 1 


order (3 2 


order (3 3 


1 x 1 


-2.0000(2) 


-1.2205(6) 


-2.9572(34) 


1x2 


-3.4488(4) 


-0.1380(13) 


-2.2088(59) 


1x3 


-4.8085(6) 


2.7795(23) 


-0.5943(72) 


1x4 


-6.1503(8) 


7.4823(38) 


-0.5589(96) 


1x5 


-7.4872(10) 


13.961(6) 


-4.460(16) 


1x6 


-8.8225(13) 


22.217(9) 


-14.684(27) 


1 x 7 


-10.157(2) 


32.247(14) 


-33.572(51) 


1x8 


-11.491(2) 


44.054(19) 


-63.501(81) 


1x9 


-12.825(2) 


57.637(26) 


-106.84(12) 


1 x 10 


-14.159(3) 


73.001(34) 


-166.05(19) 


1 x 11 


-15.492(3) 


90.135(44) 


-243.38(27) 


1 x 12 


-16.825(4) 


109.05(6) 


-341.19(39) 


1 x 13 


-18.158(4) 


129.73(7) 


-461.74(53) 


1 x 14 


-19.491(5) 


152.19(8) 


-607.52(71) 


1 x 15 


-20.825(5) 


176.42(10) 


-780.90(97) 


1 x 16 


-22.158(6) 


202.44(12) 


-984.2(1.3) 


2x2 


-5.4770(8) 


4.3363(32) 


0.0394(94) 


2x3 


-7.2455(11) 


11.630(6) 


-1.226(14) 


2x4 


-8.9557(15) 


21.711(10) 


-10.880(27) 


2x5 


-10.649(2) 


34.596(15) 


-33.645(55) 


2x6 


-12.335(2) 


50.290(22) 


-74.23(10) 


2x7 


-14.019(3) 


68.790(31) 


-137.28(16) 


2x8 


-15.701(3) 


90.107(43) 


-227.64(26) 


2x9 


-17.382(4) 


114.24(6) 


-349.97(38) 


2 x 10 


-19.062(4) 


141.19(7) 


-509.09(55) 


2 x 11 


-20.742(5) 


170.95(9) 


-709.63(76) 


2 x 12 


-22.421(6) 


203.51(12) 


-956.07(1.04) 


2 x 13 


-24.100(6) 


238.89(14) 


-1253.4(1.4) 


2 x 14 


-25.779(7) 


277.08(17) 


-1606.2(1.9) 


2 x 15 


-27.457(8) 


318.09(21) 


-2019.4(2.4) 


2 x 16 


-29.136(9) 


361.93(25) 


-2497.8(3.1) 


3x3 


-9.2216(18) 


23.225(12) 


-12.261(33) 


3x4 


-11.088(2) 


37.851(19) 


-39.256(68) 


3x5 


-12.919(3) 


55.637(29) 


-87.89(12) 


3x6 


-14.737(3) 


76.632(40) 


-164.05(19) 


3x7 


-16.549(4) 


100.85(6) 


-273.37(30) 


3x8 


-18.358(4) 


128.30(7) 


-421.87(46) 


3x9 


-20.164(5) 


158.99(9) 


-615.46(67) 


3 x 10 


-21.969(6) 


192.92(12) 


-859.71(95) 


3 x 11 


-23.773(7) 


230.09(14) 


-1160.6(1.3) 


3 x 12 


-25.576(8) 


270.49(18) 


-1524.0(1.8) 


3 x 13 


-27.379(8) 


314.12(21) 


-1955.6(2.3) 



Loop dimensions 


order p 


2 

order p 


a— 3 

order p 


3 x 14 


-29.181(9) 


361.00(25) 


-2461.2(2.9) 


3 x 15 


-30.984(10) 


411.12(30) 


-3046.9(3.7) 


3 x 16 


-32.785(12) 


464.49(35) 


-3718.3(4.7) 


4x4 


-13.051(3) 


56.867(31) 


-91.10(13) 


4x5 


-14.960(4) 


^7/~\ -1-1-1 / ^ ^ \ 

79.111(45) 


-172.59(22) 


4x6 


-16.848(4) 


104.72(6) 


c\r\r\ Oi~\/or"\ 

-290.32(35) 


4x7 


-18.725(5) 


133.74(8) 


-450.58(54) 


4x8 


-20.597(6) 


166.20(10) 


-659.90(76) 


4x9 


-22.465(7) 


202.09(13) 


-924.5(1.1) 


4 x 10 


-24.330(9) 


241.41(16) 


-1250.7(1.5) 


4 x 11 


-26.194(10) 


284.18(19) 


-1645.3(2.0) 


4 x 12 


-28.058(11) 


330.41(23) 


/ -\ -t ~t A / ' A /-"a \ 

-2114.5(2.6) 


4 x 13 


-29.918(11) 


380.09(28) 


-2665.0(3.3) 


4 x 14 


-31.779(12) 


433.19(32) 


-3302.6(4.2) 


4 x 15 


-33.640(13) 


489.77(38) 


-4034.5(5.2) 


4 x 16 


-35.501(13) 


549.80(44) 


-4866.5(6.5) 


5x5 


-16.926(5) 


105.72(7) 


-294.53(39) 


5x6 


-18.860(5) 


135.72(9) 


-460.82(57) 


5x7 


-20.778(6) 


169.20(11) 


-678.29(83) 


5x8 


-22.690(7) 


206.24(13) 


-954.0(1.1) 


5x9 


-24.596(8) 


/ -\ i /-a -| / -1 /-a \ 

246.81(16) 


-1294.0(1.6) 


5 x 10 


-26.498(9) 


290.93(20) 


-1705.6(2.1) 


5x11 


-28.398(9) 


i -\ i-x /-a -1 / 1 \ 

338.61(24) 


-2195.7(2.7) 


5 x 12 


-30.296(11) 


389.88(28) 


/ -\ ^-7 ^-7 -| -| / t-\ \ 

-2771.1(3.5) 


5 x 13 


-32.192(12) 


AAA -1 / t~\ *~\ \ 

444.71(33) 


-3438.0(4.4) 


5 x 14 


-34.088(13) 


503.06(38) 


-4202.8(5.4) 


5 x 15 


-35.982(14) 


565.02(44) 


-5073.8(6.7) 


5 x 16 


-37.878(14) 


630.53(51) 


-6056.5(8.3) 


6x6 


-20.830(6) 


170.06(11) 


-683.57(86) 


6x7 


-22.781(7) 


207.92(14) 


S~A ^-7 / -1 r~\ \ 

-965.7(1.2) 


6x8 


-24.722(8) 


249.38(17) 


-1314.9(1.5) 


6x9 


-26.654(9) 


/--\ a a -i i c~\ r\ \ 

294.41(20) 


-1737.2(2.1) 


6 x 10 


-28.582(10) 


343.08(25) 


-2240.2(2.9) 


6x11 


-30.506(10) 


395.38(29) 


-2831.0(3.6) 


6 x 12 


-32.429(12) 


451.33(34) 


t -v n^- -| /-a ^7 / j A \ 

-3516.7(4.4) 


6 x 13 


-34.350(13) 


510.92(39) 


-4303.8(5.4) 


6 x 14 


-36.269(14) 


574.13(46) 


-5198.8(6.8) 


6 x 15 


-38.188(16) 


641.05(52) 


-6211.6(8.3) 


6 x 16 


-40.104(16) 


711.56(60) 


-7345(10) 


7x7 


-24.759(8) 


250.11(18) 


-1320.6(1.7) 


7x8 


-26.722(9) 


295.93(21) 


-1751.1(2.1) 


7x9 


-28.674(10) 


345.32(25) 


-2262.4(2.9) 


7 x 10 


-30.622(11) 


398.42(31) 


-2863.8(3.8) 



Loop dimensions 


order p 


2 

order p 


order p 


7x11 


-32.566(12) 


455.16(36) 


-3561.9(4.7) 


7 x 12 


-34.504(13) 


515.62(41) 


-4364.0(5.8) 


7 x 13 


f\ /-a A A c'~\ / -t A \ 

-36.442(14) 


579.76(47) 


^ ' "\ ^-7 / ^-7 \ 

-5276.8(7.0) 


7 x 14 


-38.379(15) 


647.56(54) 


-6306.9(8.7) 


7 x 15 


-40.316(16) 


^7 -I /~\ -I -I / /"i -1 \ 

719.11(61) 


-7464(10) 


7 x 16 


-42.247(18) 


794.28(71) 


-8751(12) 


8x8 


-28.706(9) 


346.09(25) 


-2270.8(2.7) 


8x9 


-30.676(10) 


399.80(29) 


-2879.5(3.6) 


8 x 10 


-32.640(12) 


457.27(35) 


\ ^-7 / A ^-7 \ 

-3587.5(4.7) 


8x11 


-34.596(12) 


518.38(40) 


-4400.3(5.7) 


8 x 12 


-36.551(14) 


583.23(46) 


-5325.6(6.9) 


8 x 13 


-38.503(15) 


651.77(53) 


-6370.2(8.4) 


8 x 14 


-40.451(16) 


724.04(60) 


-7542(10) 


8 x 15 


A ' "\ ^ /~a -1 / -t /~a \ 

-42.401(18) 


800.10(68) 


i*~a p~ r\ i -I ' "\ \ 

-8850(12) 


8 x 16 


-44.346(19) 


879.77(77) 


-10296(14) 


9x9 


-32.660(12) 


A V t—7 l~l i~\ / O f \ 

457.82(35) 


-3593.5(4.8) 


9 x 10 


-34.639(13) 


519.60(42) 


-4415.1(6.1) 


9x11 


l"\ /-a /-a y-\ y-> / a \ 

-36.608(14) 


584.99(48) 


-5349.0(7.4) 


9 x 12 


-38.573(15) 


654.15(54) 


-6404.9(8.8) 


9 x 13 


-40.536(17) 


727.01(62) 


-7588(H) 


9 x 14 


-42.493(17) 


803.63(69) 


i*~a y^\ / -I i \ \ 

-8908(13) 


9 x 15 


-44.453(19) 


884.06(78) 


-10372(15) 


9 x 16 


-46.406(20) 


968.11(88) 


-11983(18) 


10 x 10 


-36.629(15) 


585.69(51) 


-5358.9(7.9) 


10 x 11 


-38.609(16) 


655.37(57) 


-6422.4(9.4) 


10 x 12 


-40.583(17) 


728.79(65) 


-7615(11) 


10 x 13 


-42.553(18) 


805.91(73) 


-8944(13) 


10 x 14 


-44.523(20) 


886.82(81) 


-I /~a J -1 / -1 /^a \ 

-10418(16) 


10 x 15 


-46.490(21) 


971.52(91) 


-12045(18) 


10 x 16 


-48.451(22) 


1059.8(1.0) 


-13826(22) 


11 x 11 


-40.595(18) 


729.32(66) 


^-7 /""a /™\ J / 1 1 \ 

-7624(H) 


11 x 12 


-42.581(18) 


806.97(74) 


-8964(13) 


11 x 13 


-44.557(20) 


888.30(84) 


-10445(16) 


11 x 14 


-46.530(21) 


973.41(93) 


-12080(19) 


11 x 15 


-48.504(22) 


1062.3(1.0) 


-13878(22) 


11 x 16 


-50.471(25) 


1154.8(1.2) 


-15837(26) 


12 x 12 


-44.570(21) 


888.82(85) 


-l /~\ A p~ p~ / -1 /^a \ 

-10455(16) 


12 x 13 


-46.554(23) 


/~\ f—7 A O A / s~\ r" \ 

974.34(95) 


-I r\r\r\ / -i r\ \ 

-12097(19) 


12 x 14 


-48.533(23) 


1063.6(1.0) 


-13901(22) 


12 x 15 


-50.512(24) 


1156.7(1.1) 


-15875(25) 


12 x 16 


-52.482(26) 


1253.3(1.3) 


-18017(30) 


13 x 13 


-48.541(24) 


1064.0(1.1) 


-13907(22) 


13 x 14 


-50.526(25) 


1157.4(1.2) 


-15886(25) 


13 x 15 


-52.506(26) 


1254.6(1.3) 


-18042(29) 


13 x 16 


-54.481(28) 


1355.3(1.4) 


-20375(34) 



T 

Loop dimensions 


1 

order p 


a— 2 

order p 


order p 


14 x 14 


-52.513(26) 


1254.9(1.3) 


-18050(29) 


14 x 15 


-54.496(27) 


1356.2(1.4) 


-20397(33) 


14 x 16 


-56.472(29) 


1460.9(1.6) 


-22926(39) 


15 x 15 


-56.484(28) 


1461.4(1.5) 


-22939(38) 


15 x 16 


-58.459(31) 


1570.1(1.7) 


-25672(44) 


16 x 16 


-60.433(33) 


1682.6(1.9) 


-28614(51) 



Appendix B. Unquenched n f = 2 Wilson Loops 



U(J(J\J LlllllCllolWllO 
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UI Llcl p 


1 X 1 


-Z.UUU1(1 J 
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O A 1 Q/1 ( Ci\ 

-z.41z4(^yj 


1 XX O 

i x z 


-o.44y i{Z ) 


nil op. ( i\ 
v.L16o{( ) 


1 1^QO/1 

-1.4ooz(lo J 


1 xx Q 

1 X o 


A QCidCif A\ 

-4.oUyU(4 ) 


o.l4oU(l0 ) 


U.Uoy4(^ZO ) 


1X4 


-o.iouy^o ) 


n ntxCxfi('Qn N \ 


n QQOI /K1 "\ 

-U.oozl^ol ) 


1 X/ K 

1X0 


1 A QQC\{1\ 

- 1 .4ooU( ( ) 


1/1 £ A FL f Z.\ 

14.040(0 J 


A 0/10/1 Q\ 

-4.y4y(io ) 


1 X/ £ 

1X0 


-o.ozoo^y ) 


zz.yiU^o J 


1 £ 1 £Q/QC\ 

-lo.lbo^zoj 


1 XX V 

1 X / 


-lU.loo(l ) 


QQ H/1 7/1 OA 

oo. U4 1 yi-Z J 


-oo.ooz(4o J 


1 XX Q 

1 X O 


1 1 /lOO/l "\ 

-ii.4y/(i ) 


/I /I OKO/I ft"! 
44.yoz^lO J 


-0 / .ol0(^ to) 


1 xx n 

i x y 


-LZ.oZbyZ) 


cq ^cn/oi \ 
Do.OOU^zlJ 


1 1 Q 00/1 1 \ 

-llo.UU(ll) 


1 xx 1 n 
1 X 1U 


1/11 CO/QA 

-14. loy^zj 


7/1 1 1 Q ( OQ\ 

(4:.lLo{Zo) 


1 7 /I O/l /1 7^ 

-1 /4.z4(l / J 


1 XX 1 1 

1 X 11 


i /i no f o\ 
-loA\)6{Z J 


yi. 000(^00 J 


-zoo.oo(^z4J 


1 XX 1 O 


-lo.ozo^o ) 


11U.O i [4:) 


-o04.Z0(o4 ) 


1 XX 1 Q 

1 X la 


1 o 1 c.n/q\ 


1Q1 1fi/'Cx N \ 

lol. lo^oj 


A 11 fifi/'/l 1\ 

-4 / /. oo(^4 I ) 


1 XX 1 /l 

1 X 14 


-iy.4yz^o j 


LOO. < 0{< ) 


P.OI C\P.( P.A \ 
-VZ 1 .U0^04 ) 


1 x 15 


-20.825(4) 


178.08(8) 


-804.16(85) 


1 x 16 


-22.159(4) 


204.21(10) 


-1011.6(1.1) 


2x2 


-5.4773(5) 


4.7884(23) 


0.6146(51) 


2x3 


-7.2455(8) 


12.253(5) 


-1.4455(98) 


2x4 


-8.9549(10) 


22.494(8) 


-12.438(19) 


2x5 


-10.647(1) 


35.533(12) 


-37.065(37) 


2x6 


-12.333(2) 


51.376(18) 


-80.018(69) 


2x7 


-14.016(2) 


70.022(25) 


-145.99(12) 


2x8 


-15.697(2) 


91.487(34) 


-239.76(19) 


2x9 


-17.377(3) 


115.76(5) 


-365.94(28) 


2 x 10 


-19.057(3) 


142.84(6) 


-529.34(42) 


2 x 11 


-20.736(4) 


172.74(7) 


-734.62(59) 


2 x 12 


-22.414(4) 


205.44(9) 


-986.34(81) 


2 x 13 


-24.092(5) 


240.95(11) 


-1289.4(1.1) 


2 x 14 


-25.771(5) 


279.28(14) 


-1648.7(1.5) 


2 x 15 


-27.449(6) 


320.42(16) 


-2068.6(1.9) 


2 x 16 


-29.127(7) 


364.38(19) 


-2554.1(2.4) 



Loop dimensions 


order p 


2 

order p 


a— 3 

order p 


3x3 


-9.2206(11) 


24.053(9) 


-14.087(26) 


3x4 


-ll. 085(1) 


38.856(14) 


-43.287(49) 


3x5 


-12.915(2) 


56.814(20) 


-94.774(87) 


3x6 


-14.732(2) 


77.967(29) 


A t—l A C\ A / -A A \ 

-174.34(14) 


3x7 


-16.542(3) 


102.33(4) 


-287.74(22) 


3x8 


-18.350(3) 


129.94(5) 


-440.99(33) 


3x9 


-20.154(4) 


160.77(6) 


-639.63(47) 


3 x 10 


-21.958(4) 


194.83(8) 


-889.56(66) 


3x11 


-23.760(5) 


232.12(10) 


-1196.7(9) 


3 x 12 


-25.562(5) 


272.63(12) 


-1566.5(1.2) 


3 x 13 


-27.363(7) 


316.37(15) 


-2005.1(1.6) 


3 x 14 


-29.164(7) 


363.37(18) 


-2518.7(2.1) 


3 x 15 


-30.966(7) 


A -a 6~\ r\ f C'~\ -A \ 

413.60(21) 


-3112.4(2.7) 


3 x 16 


-32.765(8) 


467.06(25) 


-3792.2(3.4) 


4x4 


-13.046(2) 


58.048(21) 


-98.171(93) 


4x5 


-14.954(2) 


80.475(30) 


-183.47(15) 


4x6 


-16.840(3) 


-a r\ t~~* r\ a / a \ 

106.24(4) 


-305.51(24) 


4x7 


-18.715(4) 


135.40(5) 


-470.56(36) 


4x8 


-20.585(5) 


168.00(7) 


-685.39(52) 


4x9 


-22.452(5) 


204.01(8) 


-956.02(71) 


4 x 10 


-24.315(5) 


243.45(10) 


-1288.9(9) 


4x11 


-26.177(6) 


286.33(12) 


-1690.7(1.2) 


4 x 12 


-28.038(7) 


332.66(14) 


-2167.6(1.6) 


4 x 13 


-29.898(7) 


382.40(17) 


-2725.9(2.1) 


4 x 14 


-31.758(7) 


435.63(20) 


-3372.5(2.7) 


4 x 15 


-33.617(8) 


492.31(24) 


AAA i ^ ^ / i ^ ^ \ 

-4113.5(3.4) 


4 x 16 


-35.477(8) 


552.42(28) 


-4955.1(4.2) 


5x5 


-16.918(3) 


107.28(4) 


-310.20(27) 


5x6 


-a c~\ nr 1 /n\ 

-18. 851(3) 


137.46(6) 


-481.94(40) 


5x7 


-20.768(4) 


-i ^ -1 -a -a / r-^\ 

171.11(7) 


-705.30(56) 


5x8 


-22.678(4) 


208.30(9) 


-987.54(80) 


5x9 


-24.582(6) 


248.99(11) 


-1334.7(1.1) 


5 x 10 


-26.482(6) 


293.22(14) 


-1753.9(1.4) 


5x11 


-28.379(7) 


341.01(16) 


-2252.0(1.8) 


5 x 12 


-30.276(8) 


392.36(19) 


-2835.7(2.4) 


5 x 13 


-32.170(8) 


447.25(23) 


-3511.6(3.1) 


5 x 14 


-34.065(9) 


505.76(26) 


-4287.1(3.8) 


5 x 15 


-35.959(10) 


567.83(31) 


-5168.9(4.7) 


5 x 16 


-37.853(12) 


633.47(35) 


-6162.9(5.8) 


6x6 


-20.820(5) 


"i r\r\ ( o\ 

172.00(8) 


-711.20(64) 


6x7 


-22.769(5) 


210.03(10) 


-1000.2(9) 


6x8 


-24.709(6) 


251.64(12) 


-1356.7(1.2) 


6x9 


-26.639(6) 


296.80(14) 


-1787.0(1.5) 


6 x 10 


-28.565(6) 


345.57(17) 


-2298.4(1.9) 



Loop dimensions 


order p 


2 

order p 


a— 3 

order p 


6 x 11 


-30.488(7) 


397.96(20) 


-2897.6(2.5) 


6 x 12 


-32.408(8) 


A f A AA / ao\ 

454.00(23) 


-3592.4(3.1) 


6 x 13 


6~\ A A a '"\ /-a / /"\ \ 

-34.326(9) 


513.66(27) 


-4389.8(4.0) 


6 x 14 


raj /-t /-~\ a A ^ / -I /~\ \ 

-36.243(10) 


577.00(31) 


-5296.2(4.8) 


6 x 15 


-38.161(11) 


643.97(36) 


-6318.2(5.8) 


6 x 16 


-40.077(12) 


714.60(40) 


-7464.1(7.0) 


7x7 


-24.746(5) 


252.41(12) 


-1363.2(1.2) 


7x8 


-26.707(6) 


298.37(15) 


-1801.5(1.6) 


7x9 


-28.658(7) 


347.89(18) 


-2321.9(2.1) 


7 x 10 


-30.604(7) 


401.08(21) 


-2932.2(2.7) 


7x11 


-32.545(8) 


457.91(24) 


-3639.4(3.3) 


7 x 12 


-34.481(9) 


518.42(27) 


-4450.8(4.1) 


7 x 13 


-36.419(10) 


582.64(32) 


-5374.8(5.0) 


7 x 14 


-38.353(10) 


650.58(36) 


-6417.0(6.0) 


7 x 15 


-40.288(H) 


722.20(41) 


n v r\ a /~* / 1—7 o \ 

-7584.6(7.3) 


7 x 16 


-42.219(11) 


797.55(46) 


AAA/ 1 f—r f C\ /"i \ 

-8886.7(8.6) 


8x8 


-28.689(7) 


348.65(19) 


-2329.4(2.2) 


8x9 


-30.658(8) 


402.50(22) 


-2948.0(2.7) 


8 x 10 


-32.619(9) 


460.03(26) 


-3665.1(3.5) 


8 x 11 


(~i A V f~7 A / t~\ \ 

-34.574(9) 


521.22(29) 


A A c\ n V f A c\\ 

-4487.5(4.2) 


8 x 12 


-36.526(10) 


586.12(33) 


-5422.8(5.1) 


8 x 13 


-38.476(11) 


654.75(38) 


-6479.3(6.2) 


8 x 14 


-40.424(11) 


^T/"\^T -I A i A 

727.14(43) 


^— /--a /--a /™\ / \ / ^— — - \ 

-7662.9(7.5) 


8 x 15 


-42.370(12) 


/~a r\ A \ c'~\ / A /"\ \ 

803.25(49) 


-8981.2(9.0) 


8 x 16 


A A A ^ -| j / h t~\ \ 

-44.314(13) 


883.16(55) 


-l /~\ A A A i -I -I \ 

-10444(11) 


9x9 


-32.639(9) 


460.64(27) 


-3672.9(3.5) 


9 x 10 


-34.615(10) 


522.46(31) 


-4504.0(4.5) 


9x11 


-36.582(11) 


587.94(34) 


a A s-\ /^a / ^ A \ 

-5448.6(5.4) 


9 x 12 


-38.546(12) 


657.15(39) 


-6514.7(6.4) 


9 x 13 


-40.507(13) 


730.09(45) 


-7710.5(7.9) 


9 x 14 


-42.464(13) 


806.83(50) 


-9042.2(9.2) 


9 x 15 


-44.420(13) 


887.30(57) 


-10517(11) 


9 x 16 


-46.373(14) 


971.61(63) 


-12147(13) 


10 x 10 


l"\ /-a /-a /-v /'~v / -1 -1 \ 

-36.602(11) 


588.54(36) 


-5457.1(5.6) 


10 x 11 


-38.581(12) 


658.28(40) 


-6532.0(6.5) 


10 x 12 


A /~\ -I i -I ' "\ \ 

-40.551(12) 


A ~\ -I ^ 7 A / A \ 

731.74(45) 


^"7 *— 7 t a y^a ^| / ^-7 ^-7\ 

-7736.1(7.7) 


10 x 13 


-42.520(13) 


808.92(51) 


-9077.9(9.4) 


10 x 14 


-44.488(14) 


889.93(57) 


-10565(11) 


10 x 15 


-46.453(15) 


974.70(64) 


-i ' "a r~\ f~\ A / -i i ai \ 

-12204(13) 


10 x 16 


-48.414(15) 


1063.31(71) 


"1 yl HAT /1 T \ 

-14005(15) 


11 x 11 


-40.566(12) 


732.25(44) 


-7744.1(7.8) 


11 x 12 


-42.549(13) 


809.94(49) 


-9094.4(9.0) 


11 x 13 


-44.522(14) 


891.34(56) 


-10590(11) 


11 x 14 


-46.496(14) 


976.61(61) 


-12240(12) 
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Loop dimensions 


order jj 


order jj 


order p 


11 x 15 


A r\ a r~ / -a \ 

-48.465(15) 


1065.59(69) 


-a A r\ A r\ Z -a v \ 

-14049(15) 


11 x 16 


-50.433(16) 


1158.43(75) 


-16030(17) 


12 x 12 


-44.535(13) 


891.87(55) 


-10600(11) 


12 x 13 


-46.515(14) 


977.47(63) 


-12257(13) 


12 x 14 


-48.495(15) 


1066.92(68) 


-14076(15) 


12 x 15 


-50.471(16) 


1160.1(8) 


-16063(17) 


12 x 16 


-52.444(16) 


1257.1(8) 


-I r~i r~\ e~\ c\ / c~\r\\ 

-18232(20) 


13 x 13 


-48.501(15) 


1067.22(72) 


-14084(16) 


13 x 14 


-50.486(16) 


1160.84(78) 


-16080(18) 


13 x 15 


-52.465(17) 


1258.13(88) 


-18252(21) 


13 x 16 


-54.441(18) 


1359.29(95) 


-20614(23) 


14 x 14 


-52.477(17) 


1258.6(9) 


-18261(20) 


14 x 15 


-54.458(18) 


1360.02(96) 


-20627(24) 


14 x 16 


-56.438(19) 


1465.3(1.0) 


-23190(27) 


15 x 15 


-56.444(20) 


1465.5(1.1) 


-23191(28) 


15 x 16 


-58.426(21) 


1574.8(1.2) 


-25960(31) 


16 x 16 


-60.408(22) 


1688.2(1.3) 


-28944(35) 
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